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We study, using both theory and molecular dynamics simulations, the relaxation dynamics of a 
microcanonical two dimensional self-gravitating system. After a sufficiently large time, a gravita- 
tional cluster of N particles relaxes to the Maxwell-Boltzmann distribution. The time to reach the 
thermodynamic equilibrium, however, scales with the number of particles. In the thermodynamic 
limit, N — ¥ oo at fixed total mass, equilibrium state is never reached and the system becomes 
trapped in a non-ergodic stationary state. An analytical theory is presented which allows us to 
quantitatively described this final stationary state, without any adjustable parameters. 



I. INTRODUCTION 

Systems interacting through long-range forces behave very differently from those in which particles interact through 
short-range potentials. For systems with short-range forces, for arbitrary initial condition, the final stationary state 
corresponds to the thermodynamic equilibrium and can be described equivalently by either microcanonical, canonical, 
or grand-canonical ensembles. On the other hand, for systems with unscreened long-range interactions, equivalence 
between ensembles breaks down [J Q . Often these systems are characterized by a negative specific heat [1|-{5| in the 
microcanonical ensemble and a broken ergodicity [H, 0]- I n the infinite particle limit, N — > oo, these systems never 
reach the thermodynamic equilibrium and become trapped in a stationary out of equilibrium state (SS) H,©]. Unlike 
normal thermodynamic equilibrium, the SS does not have Maxwell-Boltzmann velocity distribution. For finite N, 
relaxation to equilibrium proceeds in two steps. First, the system relaxes to a quasi-stationary state (qSS), in which 
it stays for time t x (N), after which it crosses over to the normal thermodynamic equilibrium with the Maxwell- 
Boltzmann (MB) velocity distribution (Iol |. In the limit N — ¥ oo, the life time of qSS diverges, t x — ¥ oo, and the 
thermodynamic equilibrium is never reached. 

Unlike the equilibrium state, which only depends on the global invariants such as the total energy and momen- 
tum and is independent of the specifics of the initial particle distribution, the SS explicitly depends on the initial 
condition. This is the case for self-gravitating systems [ll[ , confined one component plasmas [12|, EH ; geophysical 
systems [3], vortex dynamics (TB4l7|. etc [3], for which the SS state often has a peculiar core-halo structure In 
the thermodynamic limit, none of these systems can be described by the usual equilibrium statistical mechanics, and 
new methods must be developed. 

In this paper we will restrict our attention to self-gravitating systems. Unfortunately, it is very hard to study these 
systems in 3d [l9l l20j . The reason for this is that the 3d Newton potential is not confining. Some particles can gain 
enough energy to completely escape from the gravitational cluster, going all the way to infinity. In the thermodynamic 
limit, one must then consider three distinct populations: particles which will relax to form the central core, particles 
which will form the halo, and particles which will completely evaporate. Existence of three distinct classes of particles, 
makes the study of 3d systems particularly difficult. On the other hand, the interaction potential in 2d is logarithmic, 
so that all the particles remain gravitationally bound. Similar to magnetically confined plasmas the stationary state 
of a 2d gravitational system should, therefore, have a core-halo structure [Hj]. We thus expect that the insights gained 
from the study of confined plasmas might prove to be useful to understand the 2d gravitational systems. 



X 



II. THE MODEL 



Our system consists of N particles with the total mass M in a two dimensional space. At t = the particles 
are distributed over the phase space with the initial distribution fo(r,v), and then allowed to relax. Our goal is 
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to calculate the one particle distribution function /(r, v), once the relaxation process has been completed and the 
stationary state has been established. For now, we will restrict our attention to the azimuthally symmetric systems. 
The mean gravitational potential at r at time t satisfies the Poisson equation 

V 2 V = 47rGmn(r;t) (1) 

where m — M/N, n(r;t) — N j /(r, v; t)d 2 v is the particle number density, and G is the gravitational constant. It 
is convenient to define dimensionless variables by scaling lengths, velocities, potential, and energy to Lq (arbitrary 
length scale), Vq = y/2GM, V'o = 2GM and Eq = MV Q 2 = 2GM 2 , respectively. In 3d space, our system corresponds 
to infinitely long parallel rods of line density m interacting through a pair potential 4>(r) = 2Gm 2 ln(r). 



III. THE THERMODYNAMIC EQUILIBRIUM: FINITE N 



If the system has a finite number of particles, after sufficiently large time r x (iV), it will relax to thermodynamics 
equilibrium with the MB distribution function, given exactly by 



MB 



C , e -/3(v 2 /2+^(r)) 



(2) 



where C is a normalization constant, /3 — 1/T is the Lagrange multiplier used to conserve the total energy, and w(r) is 
the potential of mean force [2l|. For a gravitational system of mass M, the correlations between the particles vanish 
as N becomes large, so that cj(r) ps ip(r). Substituting Eq. ([2]) in Eq. ([T]), we obtain the classical Poisson-Boltzmann 
equation in its adimensional form 



4n 2 C 



The solution of this equation is [2S 



For large r this potential must grow as 



lim [ip(r) — ln(r)] = . 



(3) 



(4) 



(5) 



which requires that /3 = 4 and A 2 = ir 2 C/2. With these values, the distribution function Eq. ([2]) automatically 
satisfies the constraint 



J d 2 rd 2 v/(r,v) = 1, 
while the value of A is obtained from the conservation of energy 

/(r,v) = £q , 



j d 2 rd 2 v 


'v 2 ip(r)~ 




.2 2 



(6) 



(7) 



where £ o is the renormalized initial energy, see Appendix [X] In this paper we will restrict our attention to initial 
distributions of the water-bag form, 



/o(r, v) = r/Q(r m - r)Q(v m - v) 



(8) 



where Q(x) is the Heaviside step function and r\ = 1/tt r^v^. For simplicity, from now on we will measure all lengths 
in units of r m , so that r m = 1 . The renormalized energy, Appendix [S] then reduces to 



v 2 1 



(9) 



Performing the integral in Eq. ([7]) with given by Eq. (j4]), we obtain 



A 2 _ e 2(2£ -l) 



(10) 
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This provides a complete solution to the equilibrium thermodynamics of 2d self-gravitating system in a large (but 
finite) N limit. We next compare the analytical solution presented above with the full iV-body molecular dynamics 
simulation 19]. To do this we calculate the number density of particles between [r, r + dr] 

/2NX 2 r 
d 2 v/MB(r,v)= (A2+r2)2 (11) 

and the number density of particles with velocity between [v, v + dv], 

N(v) = 2irNv [ d 2 r/ A fij(r,v) = ANve' 2v2 . (12) 




FIG. 1: (a) Position and (b) velocity distributions for a system with So = —0.0433673. Solid line is the theoretical prediction 
obtained using the MB distribution function, Eq. ([2]), and the points are the result of molecular dynamics simulation with 
N = 10000 particles. 

Fig. [T]shows an excellent agreement between the theory and the simulations. It is important to stress, however, that 
to reach the MB equilibrium distribution, has required a week of CPU time, (a million dynamical times for N = 10000 
particles, see Section [VII[) . Up to the crossover time r x (N), the system remained trapped in a quasi-stationary state, 
with the one particle distribution very different from the equilibrium one. We now turn to the discussion of this 
non-equilibrium quasi-stationary state. 



IV. THE THERMODYNAMIC LIMIT 



For systems interacting through short-range potentials, the final stationary state corresponds to the thermodynamic 
equilibrium and is exactly described by the MB distribution. In spite of a popular believe that in the thermodynamic 
limit for systems with long-range unscreened interactions the mean-field becomes exact, this is not quite true. Or 
rather, this is true mathematically, but is irrelevant for real physical systems, since when N — > oo it takes an infinite 
time for such a system to relax to the thermodynamic equilibrium. What is correct, is that in the thermodynamic limit 
the dynamical evolution of a system with long-range interactions is governed exactly by the collisionless Boltzmann 
(Vlasov) equation [23| 

^ = ^ + v-V/ + --V v / = 0, (13) 

Dt Ot 771 

where / is the one particle distribution function and F is the mean force felt by a particle at position r. The MB 
distribution, together with the Poisson equation for the mean-field potential, is a stationary solution of the Vlasov 
equation. Thus, if we start with this distribution it is guaranteed to be preserved by the Vlasov dynamics. However, 
unlike for the collisional Boltzmann equation, MB distribution is not a global attractor of the Vlasov dynamics 
- an arbitrary (non-stationary) initial distribution will not evolve to the MB distribution. Thus, the collisionless 
relaxation described by the Vlasov equation is much more complex than the collisional relaxation governed by the 
Boltzmann equation for systems with short-range interactions. The final stationary state of Vlasov dynamics will 
depends explicitly on the initial particle distribution. 

Vlasov equation has an infinite number of conserved quantities, called Casimir invariants. Any local functional 
of the distribution function is a Casimir invariant of the Vlasov dynamics. In particular, if we discretize the initial 
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distribution function into surface levels with values {%}, the hypervolume = J S(f(r,v;t) — 77j)d d rd d v of each 
level will be preserved by the Vlasov flow. The evolution of the distribution function corresponds to the process of 
filamentation and proceeds ad infinitum from large to small length scale. Thus on a fine-grain scale, the evolution 
never stops and the stationary state is never reached. However, in practice, there is always a limit to the maximum 
resolution, and only a coarse-grained distribution function is available in simulations or in experiments. It is this 
coarse-grained distribution which appears in practice as the stationary state of a collisionless relaxation dynamics. 

Numerical solution of the Vlasov equation is a very difficult task. Since 1960's there has been a tremendous 
effort to find an alternative way to predict the final stationary state without having to explicitly solve the Vlasov 
equation 0, [ill, El, EM HH IHl ■ One of the first statistical approaches was proposed by Lynden-Bell and has become 
know as the violent relaxation theory . This theory is based on the assumption that there exists an efficient phase space 
mixing during the dynamical evolution. This assumption is similar to the ergodicity of the Boltzmann-Gibbs statistical 
mechanics. For systems with short-range interactions, ergodicity — although very difficult to prove explicitly — is 
almost always found to be satisfied in practice. This, however, is not the case for the efficient mixing hypothesis for 
systems with long-range interactions. In fact it was found that for most initial conditions, the phase space mixing 
is very poor. For magnetically confined plasmas, efficient mixing was found to exist only for very special initial 
conditions, and in general these systems relax to a stationary state very different from the one predicted by the 
Lynden-Bell theory. Similarly, for 3d gravitational systems, violent relaxation theory was found to work only if the 
initial distribution satisfied the, so called, virial condition (l9ll20j|. Otherwise strong particle- density wave interactions 
broke the ergodicity and resulted in a core-halo phase separation. 



A. Violent relaxation 



We first briefly review the violent relaxation theory. The basic assumption of this theory is that during the temporal 
evolution, the system is able to efficiently explore the whole of phase space. To obtain the stationary (coarse-grained) 
distribution /(r, v), the initial distribution /o(r, v) is discretized into the p levels, and the phase space is divided into 
macrocells of volume d d r d d v, which are in turn subdivided into v microcells, each of volume h d , for a d-dimensional 
system. Since Vlasov dynamics is incompressible, Df/Dt = 0, each microcell can contain at most one discretized level 
rjj . The number density of the level j inside a macrocell at (r, v) — number of microcells occupied by the level j 
divided by v — will be denoted by Pj(r, v). Note that by construction, the total number density of all levels in a 



(a) (b) 




Macro-cell 
/>,v) 

Micro-cell 
/(r,v) 



FIG. 2: Coarsening of phase-space described by the Vlasov dynamics: (a) initial and (b) final stationary state for a distribution 
with initial phase-space density rj. In this example, d = 1, p = 1 and v = 9. 

macrocell is restricted to be 

X>(r,v)<l, (14) 

3 

see Fig. [2] Using a standard combinatorial procedure [!, EI| it is then possible to associate a coarse-grained entropy 
with the distribution of {pj}. The entropy is found to be that of a p-species lattice gas, 

S = I ^^lEft^^Hft^v^ + Cl-^ft^v^Hl-^ft^v)]! , (15) 
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with the Boltzmann constant set to one. If the initial condition is of the water-bag form, Eq. ([8]) (p = 1), the 
maximization procedure is particularly simple, yielding a Fermi-Dirac distribution, 

/>, v) = i 1P (r, v) = emr J_^ } - 1 , (16) 

where e(r, v) = ^v 2 + tp(r) is the mean particle energy, and f3 and /x are the two Lagrange multipliers required by the 
conservation of the total number of particles and the total energy Eqs. (j6|) and ((TJ). 



V. VIRIAL CASES 



For a 2d self- gravitating system the virial theorem requires that (v 2 ) — 1/2, in a stationary state Appendix IB1 If the 
initial distribution does not satisfy this condition, the system will undergo strong oscillations before relaxing into the 
final stationary state in which the virial theorem is satisfied. For a water-bag initial distribution the virial condition 
reduces to the requirement that v m = 1. For future convenience, we will define the virial number for water-bag 
distributions to be fi = l/v m , so that fx = 1, when the initial distribution satisfies the virial condition. If fj, ^ 1, the 
envelope radius, defined as r e (t) = y / 2(r 2 ), will vary with time until a stationary state is achieved. Note that with 
the above definition, r e (0) = 1, as it should be. It is possible to show that the temporal evolution of r e (t) satisfies 



i 



s 2 (t) 



(17) 



r e {t) rf(t) ' 

where e 2 = 4 [(r 2 )(r 2 ) — (r • r) 2 ]. The derivation of this equation is given in Appendix [Cl For an initial water-bag 
distribution, (r(0) -r(0)) = and (r 2 (0))(r 2 (0)) = so that, if the initial distribution satisfies the virial condition, 

V m = 1 =r~ r e = 0, and the large envelope oscillations are suppressed. As was already noted for magnetically confined 
plasmas and 3d self-gravitating systems [IH, Qjl], we expect the violent relaxation theory to work well when the 
initial distribution satisfies the virial condition and there are no macroscopic envelope oscillations. To check this, we 
compare the predictions of the theory with the full iV-particle molecular dynamics simulations. At t = 0, particles are 
distributed over the phase space in accordance with the water-bag distribution ([5]) which satisfies the virial condition, 
fj, = 1. We then numerically solve the Poisson equation (TTJ), with the distribution function given by equation (|16[) 




FIG. 3: Position (a) and velocity (b) distributions for a system satisfy the virial condition. The solid line is the theoretical 
prediction obtained using the distribution function of Eq. (|16|l , while the points are the result of molecular dynamics simulation 
with N = 10000 particles. 



and compare the results with the molecular dynamics simulations. As can be seen from Fig. [3] there is a reasonably 
good agreement between the theory and the simulations. However, if the virial condition is not met exactly, one 
notices a deviation in the tail region of the particle distribution, see inset of Fig. @] For initial distributions with 
(i significantly different from 1, there is a clear qualitative change in the SS distribution function. In this case the 
original homogeneous cluster, separates into a high density core region surrounded by a diffuse halo, Fig. [5] — the 
violent relaxation theory fails completely and a new approach must be developed [H, [ljl H(| . 
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FIG. 4: (a) Position distribution for a nearly virial self-gravitating system with initial /i = 1.2. Solid line is the theoretical 
prediction obtained using the distribution function of Eq. H16|) . while the points are the result of molecular dynamics simulation 
with N = 10000 particles. 




FIG. 5: Position (a) and velocity (b) distributions for a system with fi = 1.7. Solid line is the prediction of the violent relaxation 
theory Eq. (|16|) . while points are the result of molecular dynamics simulation with N — 10000 particles. 



VI. CORE-HALO DISTRIBUTION 



The failure of the violent relaxation theory is a consequence of the inapplicability of the efficient mixing hypothesis 
to strongly oscillating gravitational systems, Figs. 0] and [5j Density oscillations excite parametric resonances which 
favor some particles to gain a lot of energy at the expense of the rest. The resulting particle- wave interactions are 
a form of non-linear Landau damping which allows some particles to escape from the main cluster to form a diffuse 
halo. The process of evaporation will continue as long as the oscillations of the core persist. Oscillations will only stop 
when the core exhaust all of its free energy, and its effective temperature drops to T « 0, — > oo in Eq. (flU)) . Note 
that because of the incompressibility restriction imposed by the Vlasov dynamics Eq. (|T4")) . the core can not freeze - 
collapse to the minimum of the potential energy. Instead, the distribution function of the core particles progressively 
approaches that of a fully degenerate Fermi gas [l9| , 

/core(r,v) = 7ye(e F -e(r,v)) (18) 

where tp is the effective Fermi energy. The final stationary state of the cluster will then correspond to a cold core 
surrounded by a high energy diffuse halo, 

/(r,v) = r;e( eF -e(r,v))+xe( e (r,v)-e F )e( efl -e(r,v)), (19) 

where en is the energy of the one particle resonance. The parameter \ an d the effective Fermi energy tp are determined 
using the conservation of particle number and energy. The extent and the location of the parametric resonance can 
be calculated using the canonical perturbation theory j27j. In Fig. [5] we show the Poincare section of a test particle i 
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moving under the action of an oscillating potential calculated using the envelope equation (fl7|) . 

" for n {t) < r e (t) 



hit) 



L 



'(*) 



~7k) for n(t)>r e (t) 



with, 



Te{t) rj(t) ' 



(20) 



(21) 



where L t = |r< x Vj| is the modulus of the test particle angular momentum, conserved by the dynamics, and e(t) is 
fixed at its initial value e(t) = e = v m = 1//J. The resonant orbit is the outermost curve of the Poincare plot Fig. El 
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FIG. 6: (a) Poincare plot of a test particle in a oscillating potential, Eq. (|20l) with 10 different initial conditions, plotted when 
the envelope is at its minimum, (b) N-particle simulation for a non-virial system with /i = 1.2. An excellent agreement is found 
between the extent of the halo in N-particle simulation and the one particle resonant orbit shown in the Poincare plot. 

The first resonant particles move in an almost simple harmonic motion with energy e R — ln(r^), where r R is the 
intersection of the resonant trajectory with the v = axis. 

Empirically we find that the location of the one particle resonance r R for values of |/x — 1| > 0.1 is very well 
approximated by a simple expression [28[ 

r« = 2 (1 + 1 IhmD/m ■ (22) 

As the relaxation proceeds, the oscillating core becomes progressively colder, while a halo of highly energetic particles is 
formed. As more and more particles are ejected from the core, their motion becomes chaotic, and a halo distribution 
becomes smeared out. Similar to what happens for magnetically confined plasmas, we find that the distribution 
function of a completely relaxed halo is very well approximated by the Heaviside step function 0(e# — e(r, v)). For 
notational simplicity, from now on, we will drop the over-bar on the distribution function / (r, v), but it should always 
be kept in mind that / is stationary only within the coarse graining procedure described above. 



VII. ANALYTICAL SOLUTION TO THE CORE-HALO PROBLEM 

In order to obtain the density and the velocity distribution after the SS state is achieved, we solve the Poisson 
equation JT]) 

V 2 V(r) = 27ry"/(r,v)d 2 v, (23) 

with the constraints ([6]) and (JT]). Since the initial mass distribution has the azimuthal symmetry, the potential must 
have the form 

^(r) = ipcore(r)Q(r c -r)+ip ha i (r)Q(r-r c )Q(r R -r)+tp out {r)Q(r -r R ). (24) 

Substituting this into Poisson equation and noting that cf — ip{ r c) we obtain 

Vw(r) = e R + C 1 [(T)/x-l)MK) + Mr*)} (25) 
i>haio(r) = e R + C 2 J (r**) + C 3 Y (r**) (26) 
VwM = ln(r) , (27) 
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where Jo and lo are the Bessel functions of the first type and of order 0; r c is the core radius; r* = 2nr^fr\ and 
r** = l-nr^fx- The integration constants Ci.2,3 and the value of r c can be determined using the continuity of the 
potential and of the gravitational field. The parameter x can then be obtained using the conservation of energy. Once 
Ci.2,3 are calculated, see Appendix iDl we are left with just two equations for r c and x, 

£(r c ,X)-£o = (28) 

*P'core( r c) = f ha l (rc) , 

where prime denotes the derivative with respect to r and 

Ft x _ **x(v - X)J2 (r*) rj [Y (r**) Jo (rg) - Jo (r e ** ) Yp (r%)f 

[rc,X) ~ 2 47?J (r*) ■ [ZJ) 

(30) 

This completely determines the distribution function of the final stationary state achieved by a self-gravitating 
system when its initial distribution deviates from the virial condition. In Fig. [71 we compare the predictions of the 
theory with the molecular dynamics simulations. An excellent agreement is found without any adjustable parameters. 




FIG. 7: Position (a) and velocity (b) distributions for a system with fx — 1.7. Solid line is the theoretical prediction obtained 
using the distribution function, Eq. (|19[l . and the points are the result of molecular dynamics simulation with TV = 10000 
particles. 

Finally, we explore the life time t x (N) of a qSS of a self-gravitating system with a finite number of particles. To 
do this we define the crossover parameter 



1 f c 



[N(r,t)-N ch (r)] 2 dr (31) 



where N(r, t) is the number density of particles inside shells located between r and r + dr at each time of simulation 
t and N c h(r) = 2-nNr J f c h (r, v) d 2 v where f c h (r, v) is the stationary distribution given by Eq. (|T5)> . The dynamical 
time scale is defined as td = r m /y/2GM. In Fig. [8^ we plot the value of for systems with different number of 
particles. Fig. [5Jd shows that if we scale the time with r x = N 7 T£>, where 7 = 1.35, all the curves collapse onto one 
universal curve, showing the divergence of the crossover time in the thermodynamic limit. It is interesting to note 



that for a Hamiltonian Mean-Field (HMF) model, t x (N) was found to diverge with the exponent 7 = 1.7 31], while 
for a virial 3d self-gravitating system the exponent was found to be 7 ~ 1. Unfortunately, at the moment there is no 
theory which allows us to predict this exponents a priory. 
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FIG. 8: (a) for different number of particles in the system. After relaxing into the qSS, the system crosses over to MB 
distribution after time r x (N). Inset (a) shows that the relaxation to core-halo state takes approximately t ft* 2000td and does 
not depend on the number of particles in the system. When the time is scaled with r x (N) all the data in (a) (for large times) 
collapse onto one universal curve (b) . 
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VIII. CONCLUSIONS 

We have studied the thermodynamics of 2d self-gravitating system in the microcanonical ensemble. It was shown 
that the gravitational clusters containing finite number of particles relax to the equilibrium state characterized by 
the MB distribution. Prior to achieving the thermodynamic equilibrium, however, these systems become trapped in 
a quasi-stationary state, where they stay for time r x , which diverges as N 1 - 35 for large N. Thus, in the limit N —> oo 
at fixed total mass M, thermodynamic equilibrium can not be reached in a finite time. A new approach, based of the 
conservation properties of the Vlasov dynamics and on the theory of parametric resonances, is formulated and allows 
us to quantitatively predict the one particle distribution function in the non-equilibrium stationary state. Finally, it 
is curious to consider what will happen to a self-gravitating system in a contact with a thermal bath — the canonical 
ensemble. In Appendix [Bj it is shown that for a 2d self-gravitating system a stationary state is possible, if and only if, 
(v 2 )=l/2, i.e. when the kinetic temperature is T = 1/4. If such system is put into contact with a thermal bath which 
has T > 1/4 there will be a constant heat flux from the reservoir into the system. This heat will be converted into the 
gravitational potential energy — since the kinetic energy is fixed by the virial condition — making the cluster expand 
without a limit. Conversely if the bath temperature is T < 1/4, the heat flux will be from the system into the bath. 
Again, since the system can only exist in a stationary state if T = 1/4, the energy for the heat flux can come only from 
the gravitational potential. In this case the gravitational cluster will contract without a limit, concentrating all of its 
mass at the origin. Thus, in the canonical ensemble no thermodynamic equilibrium is possible, unless the reservoir 
is at exactly T = 1/4. We hope that the present work will also help shed new light on the collisionless relaxation in 
3d self-gravitating systems. Unfortunately the 3d problem is significantly more difficult, since besides the core-halo 
formation, one must also account for the particles evaporating to infinity. 

This work was partially supported by the CNPq, INCT-FCx, and by the US-AFOSR under the grant FA9550-09- 
1-0283. 
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Appendix A: The Total Energy 

The gravitational potential energy U of a system is 

[/ = -—/ (-V^) 2 d 2 r (Al) 
4?r J 

where the integration extends over all space. Unlike in 3d, the gravitational potential of a 2d system diverges at 
infinity. Therefore, some care must be taken with the limits. Performing the integration by parts we obtain 

U=\f Uir)- lim ^(r )W,v;t)d 2 rd 2 v, (A2) 



2 . 

where rg is the radius of the bounding sphere. From Eq. ([SJ, ip{ r o) = hi(ro), and the total energy is given by 

E = j ^ + ^^/(r,v;t)d 2 rd 2 v-i r hm o ln(r ) (A3) 

The divergence in the last term is common to all energy calculations in 2d. For example, if at t = the system is 
distributed with a water-bag distribution Eq. (|5J|, its energy is 

E Q = ^ - I + \ In (r m ) - \ lim ln(r ) . (A4) 

The gravitational self-energy in 2d is always divergent. However, since this divergence is always the same, it can be 
easily renormalized away. We simply add an infinite constant, \ linv^oo ln(ro), to all gravitational self-energies. The 
renormalized (finite) energy £ of a self-gravitating system is then 

^/(y + ^)/(r,v )d 2 rd 2 v (A5) 

which for a water-bag distribution Eq. © becomes 

£ = V -f-\ + \Hr m ). (A6) 



Appendix B: The Virial Theorem 

The Hamiltonian of a general self-confined system is given by 



where pi is the momentum of particle i, and V(ri — Tj) is the interaction potential. The virial function / is defined as 

1 = (52 riP ^ ■ 

i 

Taking the time derivative, and using Hamilton's equations we obtain 

,2 



1 m t— 1 dr 

i 

where, 



V 2 
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If V is an homogeneous function of order p 

V(r) = \- p V (Ar) , 

Euler's theorem requires that 

i 

For a stationary state d//dt = 0, and we obtain the usual result 2K, = p{V), where K, is the mean kinetic energy. 

In two dimensions, V is a sum of logarithms and Euler equation does not apply directly. However, we can still 
derive a 2d virial theorem by writing the inter-particle interaction potential as ^(r) = 2Gm 2 lim p _>.o(|r| p )/p, which is 
a logarithm plus an infinite constant. This is a homogeneous function of order p=0, so we can use the Euler theorem 
to write 

Gm 2 N(N- 1) = J2 Ti 7rV ( B2 ) 

i 

Note that since the right hand side of this equation only contains derivative of the potential, this expression is 
equally valid for p = and for the logarithmic potential. Substituting Eq. (|B2[) into Eq. (|B1[) , we arrive at 2d virial 
theorem [30j | . 

(« 2 > = GM^-i (B3) 



In the thermodynamic limit, and after rescaling the velocity to put everything into adimensional form, we obtain the 
result (v 2 ) — 1/2, quoted in Section IVl 
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Appendix C: The Envelope Equation 

We define the rms radius of the mass distribution as R = yj (r 2 ). Deriving twice with respect to time we obtain 



R = 



(r 2 )(r 2 



R 3 



R 



(CI) 



This reduces to 



R = 



4i? 3 



R 



(C2) 



where e 2 = 4 ((r 2 )(r 2 ) — (r • r) 2 ) is the emittance which commonly appears in plasma physics [29(. The last term can 
be simplified using the Poisson equation ([!]) in its dimensionless form, 



(r • r) = J r • r /(r, v,t)d 2 vd 2 v 



1 

2^ 



= - I r 



r • r V tjjd r 

^ 2 
— Vipdr 
or 

dip d ( dtp 



dr dr V dr 



dr 



dr 



0_ 

dr 



chp_ 
dr 



i ( d^py 

1 ™ r 9r j 



(C3) 



which can be obtained directly using Eq. ([5]), 



■1/2. 



(C4) 



For a water bag initial distribution flSJ) we define the envelope radius as r e = Ry/2, so that for t = r e (0) = 1, and 
rewrites (IC2I) as 



r e (*) 



1 



r e (t) rf(t) ' 



(C5) 



This is the envelope equation. If initially, (w 2 ) =1/2 then f e = and the envelope will not oscillate. This is precisely 
the virial condition. 
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Appendix D: Potential for Core-Halo Distribution 

Integrating the core-halo distribution function over velocity, the dimensionless Poisson equation (|23[) becomes 



rj(e F - ijj) + x(e_R - e F ) for ip < e F 
V 2 -0 = 4?r 2 { xUr ~ V) for e F < ip < e R (Dl) 
for -0 > e R . 

We define ip cor e for ty < £f, Whalo for e F < ip < e R , and ipout for W > e R- Changing variables r* = 2-Kr^fr\ and 
r** = we can rewrite (|D1[) as 

/ " i Wcore , / , X / \ 

Weave H r~ + VWe = £F + "(Cfl - e F ) 

r* T) 

^ + ^=0 (D2) 

The solution of the first two of these equations can be written in terms of the Bessel functions of first type and of 
order 0, 

Vcore{r) = e F + ^{e R -e F ) + CxJ {r*) + C v Y {r*) (D3) 
V 

VhaUr) = e R + C 2 J (r**)+C 3 Y (r**) (D4) 

The last equation is solved by 

Vout(r) = C A In r + C v , (D5) 

where {Cj} are the integration constants. The regularity of solution at the origin and Eq.([5]) require that Cy = 0, 
CV = and C4 = 1, respectively. The potential reduces to 

Veore(r) = C R + d [fa/x ~ 1) J «) + J (f*)] (D6) 

</Wr) = e^ + C 2 Jo(r**) + C 3 ro(r**) (D7) 
Vout(r) = ln(r) , (D8) 

where we have defined r c such that £p = ip(r c ). The others requirements to continuity of the potential and its 
derivative are, 

Veore{r c ) - 1phalo{r c ) = 
Vhalo(r R ) - Ipout(rR) = 

Solving these equations yields the integration constants Ci.2,3 as a function of (r c ,x) an d t ne parameters (e R ,r]), 

Ci = 7T X (y (27rr e yx) Jo (2^r fl ,^x) - J (2nr c ^x) Y (2nr R ^x)) 

2rjJ (2irr cy /rf) 

c 2 = _^pVxl (D10) 

C 3 = ^ , (Dll) 

The remaining equation of continuity of tp'(r) at r c and the conservation of energy will determine r c and x, see 
Eq.®. 
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